# TODO: Add comment
# 
# Author: yaping
###############################################################################
##parameters: working_dir, file_one, file_two, NDR_p_value_threshold, annotation_file, overlap_or_not_with_annotation_location, output_bed_file_in_different_combination_with_gene_names
##each column's content:  chr	start	end	aveMethyWind	gchNumWind	gchDepthWind	gchCTdepthWind	gchDotWind	sigValue	aveMethyWindLinker	wcgMethyWind	wcgNumWind


for (e in commandArgs(TRUE)) {
	ta = strsplit(e,"=",fixed=TRUE)
	if(! is.na(ta[[1]][2])) {
		if(ta[[1]][1] == "wd"){
			wd<-ta[[1]][2]  ## directory to make plot file
			
		}
		if(ta[[1]][1] == "file1"){
			file1<-ta[[1]][2]
		}
		if(ta[[1]][1] == "sample1"){
			sample1<-ta[[1]][2]
		}
		if(ta[[1]][1] == "file2"){
			file2<-ta[[1]][2]
		}
		if(ta[[1]][1] == "sample2"){
			sample2<-ta[[1]][2]
		}
		if(ta[[1]][1] == "pValue"){
			pValue<-ta[[1]][2]
		}
		if(ta[[1]][1] == "anno"){
			anno<-ta[[1]][2]
		}
		if(ta[[1]][1] == "annoFormat"){
			annoFormat<-ta[[1]][2]
		}
		if(ta[[1]][1] == "overlap"){
			overlap<-as.numeric(ta[[1]][2])  ##0: no overlap; 1: overlap
		}
		if(ta[[1]][1] == "outputBed"){
			outputBed<-as.numeric(ta[[1]][2])  ## 1: output; others: no
		}
		if(ta[[1]][1] == "useWcg"){
			useWcg<-as.numeric(ta[[1]][2])  ## 1: use cg methylation value; others: no
		}
		if(ta[[1]][1] == "dataPoint"){
			dataPoint<-as.numeric(ta[[1]][2])  ## 1: use cg methylation value; others: no
		}
		
	}
}
setwd(wd)
library(Ringo)

Obj1<-read.table(file1,sep="\t",header=T)
Obj2<-read.table(file2,sep="\t",header=T)

Obj1<-Obj1[!is.na(Obj1[,9]) & !is.na(Obj1[,4]) & Obj1[,8] > dataPoint,]
Obj2<-Obj2[!is.na(Obj2[,9]) & !is.na(Obj2[,4]) & Obj2[,8] > dataPoint,]

Obj1 <-cbind(Obj1, p.adjust(Obj1[,9],method="BH"))
Obj2 <-cbind(Obj2, p.adjust(Obj2[,9],method="BH"))



Obj1Data<-array()
Obj2Data<-array()

for(i in c(1:22)){
	chrName=paste("chr",i,sep="")
	tmpObj1<-Obj1[Obj1[,1] %in% chrName,]
	tmpObj2<-Obj2[Obj2[,1] %in% chrName,]
	
	Obj1LocData<-data.frame(chr= tmpObj1[,1], start= tmpObj1[,2], end= tmpObj1[,3])
	Obj2LocData<-data.frame(chr= tmpObj2[,1], start= tmpObj2[,2], end= tmpObj2[,3])
	
	Obj1Obj2tssOverlapped<-as.matrix(regionOverlap(Obj1LocData, Obj2LocData))
	
	
	Obj1Data<-rbind(Obj1Data,cbind(tmpObj1, rowSums(Obj1Obj2tssOverlapped)))
	Obj2Data<-rbind(Obj2Data,cbind(tmpObj2, rowSums(t(Obj1Obj2tssOverlapped))))
	
}
Obj1Data<-Obj1Data[2:length(Obj1Data[,1]),]
Obj2Data<-Obj2Data[2:length(Obj2Data[,1]),]

Obj1NDRInObj1<-Obj1Data[Obj1Data[,13]< pValue,]
Obj2NDRInObj2<-Obj2Data[Obj2Data[,13]< pValue,]

Obj1NDR<-Obj1Data[Obj1Data[,14]!=0 & Obj1Data[,13]< pValue,]
Obj2NDR<-Obj2Data[Obj2Data[,14]!=0 & Obj2Data[,13]< pValue,]

sharedWind<-Obj1Data[Obj1Data[,14]!=0,]

Obj1LocNDR<-data.frame(chr= Obj1NDR$chr, start= Obj1NDR$start, end= Obj1NDR$end)
Obj2LocNDR<-data.frame(chr= Obj2NDR$chr, start= Obj2NDR$start, end= Obj2NDR$end)

Obj1Obj2NDRoverlapped<-as.matrix(regionOverlap(Obj1LocNDR, Obj2LocNDR))
NumObj1Obj2NDRoverlapped<-Obj1Obj2NDRoverlapped[rowSums(Obj1Obj2NDRoverlapped)!=0,]


if(!is.null(anno)){
	
	annoFile<-read.table(anno,sep="\t",header=F)
	annoLoc<-data.frame(chr= annoFile[,1], start= annoFile[,2], end=annoFile[,3])
	if(annoFormat == "gtf"){
		annoLoc<-data.frame(chr= annoFile[,1], start= annoFile[,4], end=annoFile[,5])
	}
	callableWindowObj1Data=Obj1Data
	callableWindowObj2Data=Obj2Data
	callableWindowLocObj1<-data.frame(chr= callableWindowObj1Data[,1], start= callableWindowObj1Data[,2], end=callableWindowObj1Data[,3])
	callableWindowLocObj2<-data.frame(chr= callableWindowObj2Data[,1], start= callableWindowObj2Data[,2], end=callableWindowObj2Data[,3])
	
	
	Obj1DataAnnoInfo<-array()
	Obj2DataAnnoInfo<-array()
	
	for(i in c(1:22)){
		chrName=paste("chr",i,sep="")
		tmpObj1<-callableWindowObj1Data[callableWindowObj1Data[,1] %in% chrName,]
		tmpObj2<-callableWindowObj2Data[callableWindowObj2Data[,1] %in% chrName,]
		
		tmpAnno<-annoFile[annoFile[,1] %in% chrName,]
		
		Obj1LocData<-data.frame(chr= tmpObj1[,1], start= tmpObj1[,2], end= tmpObj1[,3])
		Obj2LocData<-data.frame(chr= tmpObj2[,1], start= tmpObj2[,2], end= tmpObj2[,3])
		
		annoLocTmp<-data.frame(chr= tmpAnno[,1], start= tmpAnno[,2], end= tmpAnno[,3])
		
		if(annoFormat == "gtf"){
			annoLocTmp<-data.frame(chr= tmpAnno[,1], start= tmpAnno[,4], end=tmpAnno[,5])
		}
		
		Obj1AnnoOverlapped<-as.matrix(regionOverlap(Obj1LocData, annoLocTmp))
		Obj2AnnoOverlapped<-as.matrix(regionOverlap(Obj2LocData, annoLocTmp))
		
		Obj1DataAnnoInfo<-rbind(Obj1DataAnnoInfo,cbind(tmpObj1, rowSums(Obj1AnnoOverlapped)))
		Obj2DataAnnoInfo<-rbind(Obj2DataAnnoInfo,cbind(tmpObj2, rowSums(Obj2AnnoOverlapped)))
		
	}
	Obj1DataAnnoInfo<-Obj1DataAnnoInfo[2:length(Obj1DataAnnoInfo[,1]),]
	Obj2DataAnnoInfo<-Obj2DataAnnoInfo[2:length(Obj2DataAnnoInfo[,1]),]
	
	Obj1Anno<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,15]!=0,]
	Obj2Anno<-Obj2DataAnnoInfo[Obj2DataAnnoInfo[,15]!=0,]
	
	Obj1NDRannoInObj1<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,15]!=0 & Obj1DataAnnoInfo[,13]<pValue,]
	Obj2NDRannoInObj2<-Obj2DataAnnoInfo[Obj2DataAnnoInfo[,15]!=0 & Obj2DataAnnoInfo[,13]<pValue,]
	
	Obj1NDRanno<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,14]!=0 & Obj1DataAnnoInfo[,15]!=0 & Obj1DataAnnoInfo[,13]<pValue,]
	Obj2NDRanno<-Obj2DataAnnoInfo[Obj2DataAnnoInfo[,14]!=0 & Obj2DataAnnoInfo[,15]!=0 & Obj2DataAnnoInfo[,13]<pValue,]
	
	sharedWindAnno<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,14]!=0 & Obj1DataAnnoInfo[,15]!=0,]
	
	if(useWcg){
		Obj1Anno<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,16]!=0 & (Obj1DataAnnoInfo[,11]<0.2 | is.na(Obj1DataAnnoInfo[,11])),]
		Obj2Anno<-Obj2DataAnnoInfo[Obj2DataAnnoInfo[,16]!=0 & (Obj2DataAnnoInfo[,11]<0.2 | is.na(Obj2DataAnnoInfo[,11])),]
		
		Obj1NDRannoInObj1<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,16]!=0 & Obj1DataAnnoInfo[,13]<pValue & (Obj1DataAnnoInfo[,11]<0.2 | is.na(Obj1DataAnnoInfo[,11])),]
		Obj2NDRannoInObj2<-Obj2DataAnnoInfo[Obj2DataAnnoInfo[,16]!=0 & Obj2DataAnnoInfo[,13]<pValue & (Obj2DataAnnoInfo[,11]<0.2 | is.na(Obj2DataAnnoInfo[,11])),]
		
		Obj1NDRanno<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,14]!=0 & Obj1DataAnnoInfo[,15]!=0 & Obj1DataAnnoInfo[,16]!=0 & Obj1DataAnnoInfo[,13]<pValue & (Obj1DataAnnoInfo[,11]<0.2 | is.na(Obj1DataAnnoInfo[,11])),]
		Obj2NDRanno<-Obj2DataAnnoInfo[Obj2DataAnnoInfo[,14]!=0 & Obj2DataAnnoInfo[,15]!=0 & Obj2DataAnnoInfo[,16]!=0 & Obj2DataAnnoInfo[,13]<pValue & (Obj2DataAnnoInfo[,11]<0.2 | is.na(Obj2DataAnnoInfo[,11])),]
		
	}
	
	if(overlap==0){
		Obj1NDRanno<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,14]!=0 & Obj1DataAnnoInfo[,15]==0 & Obj1DataAnnoInfo[,13]<pValue,]
		Obj2NDRanno<-Obj2DataAnnoInfo[Obj2DataAnnoInfo[,14]!=0 & Obj2DataAnnoInfo[,15]==0 & Obj2DataAnnoInfo[,13]<pValue,]
		sharedWindAnno<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,14]!=0 & Obj1DataAnnoInfo[,15]==0,]
		Obj1Anno<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,15]==0,]
		Obj2Anno<-Obj2DataAnnoInfo[Obj2DataAnnoInfo[,15]==0,]
		Obj1NDRannoInObj1<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,15]==0 & Obj1DataAnnoInfo[,13]<pValue,]
		Obj2NDRannoInObj2<-Obj2DataAnnoInfo[Obj2DataAnnoInfo[,15]==0 & Obj2DataAnnoInfo[,13]<pValue,]
		
		if(useWcg){
			Obj1NDRanno<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,14]!=0 & Obj1DataAnnoInfo[,15]!=0 & Obj1DataAnnoInfo[,16]==0 & Obj1DataAnnoInfo[,13]<pValue & (Obj1DataAnnoInfo[,11]>=0.2 | is.na(Obj1DataAnnoInfo[,11])),]
			Obj2NDRanno<-Obj2DataAnnoInfo[Obj2DataAnnoInfo[,14]!=0 & Obj2DataAnnoInfo[,15]!=0 & Obj2DataAnnoInfo[,16]==0 & Obj2DataAnnoInfo[,13]<pValue & (Obj1DataAnnoInfo[,11]>=0.2 | is.na(Obj1DataAnnoInfo[,11])),]
			
			Obj1Anno<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,16]==0 & (Obj1DataAnnoInfo[,11]>=0.2 | is.na(Obj1DataAnnoInfo[,11])),]
			Obj2Anno<-Obj2DataAnnoInfo[Obj2DataAnnoInfo[,16]==0 & (Obj1DataAnnoInfo[,11]>=0.2 | is.na(Obj1DataAnnoInfo[,11])),]
			
			Obj1NDRannoInObj1<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,16]==0 & Obj1DataAnnoInfo[,13]<pValue & (Obj1DataAnnoInfo[,11]>=0.2 | is.na(Obj1DataAnnoInfo[,11])),]
			Obj2NDRannoInObj2<-Obj2DataAnnoInfo[Obj2DataAnnoInfo[,16]==0 & Obj2DataAnnoInfo[,13]<pValue & (Obj1DataAnnoInfo[,11]>=0.2 | is.na(Obj1DataAnnoInfo[,11])),]
			
		}
		
	}
	
	Obj1LocNDRanno<-data.frame(chr= Obj1NDRanno$chr, start= Obj1NDRanno$start, end= Obj1NDRanno$end)
	Obj2LocNDRanno<-data.frame(chr= Obj2NDRanno$chr, start= Obj2NDRanno$start, end= Obj2NDRanno$end)
	
	Obj1Obj2NDRoverlappedAnno<-as.matrix(regionOverlap(Obj1LocNDRanno, Obj2LocNDRanno))
	NumObj1Obj2NDRoverlappedAnno<-Obj1Obj2NDRoverlappedAnno[rowSums(Obj1Obj2NDRoverlappedAnno)!=0,]
	
}


cat("Total callable window of Obj1:")
cat(dim(Obj1)[1])
cat("\n")
cat("Total callable window of Obj2:")
cat(dim(Obj2)[1])
cat("\n")
cat("Total NDR of Obj1:")
cat(dim(Obj1NDRInObj1)[1])
cat("\n")
cat("Total NDR of Obj2:")
cat(dim(Obj2NDRInObj2)[1])
cat("\n")
cat("NDR percentage in callable window of Obj1:")
cat(dim(Obj1NDRInObj1)[1]/dim(Obj1)[1])
cat("\n")
cat("NDR percentage in callable window of Obj2:")
cat(dim(Obj2NDRInObj2)[1]/dim(Obj2)[1])
cat("\n")
cat("Number of shared callable window:")
cat(dim(sharedWind)[1])
cat("\n")
cat("NDR of Obj1 in shared callable window:")
cat(dim(Obj1NDR)[1])
cat("\n")
cat("NDR of Obj2 in shared callable window:")
cat(dim(Obj2NDR)[1])
cat("\n")
cat("NDR overlapped between Obj1 and Obj2:")
cat(dim(NumObj1Obj2NDRoverlapped)[1])
cat("\n")
cat("percentage of NDR in Obj1 overlapped between Obj1 and Obj2:")
cat(dim(NumObj1Obj2NDRoverlapped)[1]/dim(Obj1NDR)[1])
cat("\n")
cat("percentage of NDR in Obj2 overlapped between Obj1 and Obj2:")
cat(dim(NumObj1Obj2NDRoverlapped)[1]/dim(Obj2NDR)[1])
cat("\n")

if(!is.null(anno)){
	cat("Total callable window of Obj1 in annotation file:")
	cat(dim(Obj1Anno)[1])
	cat("\n")
	cat("Total callable window of Obj2 in annotation file:")
	cat(dim(Obj2Anno)[1])
	cat("\n")
	cat("Total NDR of Obj1 in annotation file:")
	cat(dim(Obj1NDRannoInObj1)[1])
	cat("\n")
	cat("Total NDR of Obj2 in annotation file:")
	cat(dim(Obj2NDRannoInObj2)[1])
	cat("\n")
	cat("NDR percentage in callable window of Obj1 in annotation file:")
	cat(dim(Obj1NDRannoInObj1)[1]/dim(Obj1Anno)[1])
	cat("\n")
	cat("NDR percentage in callable window of Obj2 in annotation file:")
	cat(dim(Obj2NDRannoInObj2)[1]/dim(Obj2Anno)[1])
	cat("\n")
	cat("Number of shared callable window in annotation file:")
	cat(dim(sharedWindAnno)[1])
	cat("\n")
	cat("NDR of Obj1 in shared callable window in annotation file:")
	cat(dim(Obj1NDRanno)[1])
	cat("\n")
	cat("NDR of Obj2 in shared callable window in annotation file:")
	cat(dim(Obj2NDRanno)[1])
	cat("\n")
	cat("NDR overlapped between Obj1 and Obj2 in annotation file:")
	cat(dim(NumObj1Obj2NDRoverlappedAnno)[1])
	cat("\n")
	cat("percentage of NDR in Obj1  overlapped between Obj1 and Obj2 in annotation file:")
	cat(dim(NumObj1Obj2NDRoverlappedAnno)[1]/dim(Obj1NDRanno)[1])
	cat("\n")
	cat("percentage of NDR in Obj2  overlapped between Obj1 and Obj2 in annotation file:")
	cat(dim(NumObj1Obj2NDRoverlappedAnno)[1]/dim(Obj2NDRanno)[1])
	cat("\n")
}

if(outputBed == 1){
	#filename=paste(sample1,annoName,"pValue",pValue,"DP",dataPoint,"GCH",gchNum, "summary.txt",sep=".")
	#summary<-cbind(dim(Obj1)[1],dim(Obj1NDRInObj1)[1],dim(Obj1Anno)[1],dim(Obj1Anno)[1]/dim(Obj1)[1],dim(Obj1NDRannoInObj1)[1],dim(Obj1NDRannoInObj1)[1]/dim(Obj1Anno)[1],(dim(Obj1Anno)[1] + dim(Obj1NDRInObj1)[1])/dim(Obj1Anno)[1],(dim(Obj1Anno)[1] + dim(Obj1NDRInObj1)[1])/dim(Obj1NDRInObj1)[1])
	#colnames(summary)<-c("Obj1_callableWindow","Obj1_NDRs","Anno_in_Obj1_callableWindow","Perc_Anno_in_Obj1_callableWindow","Obj1_NDR_overlap_anno","Perc_Obj1_NDR_overlap_anno","(anno+NDR)/anno","(anno+NDR)/NDR")
	#write.table(summary,filename,quote=F,sep="\t",row.names=F,col.names=T)
	library(ChIPpeakAnno)
	data(TSS.human.GRCh37)
	Obj1NDRoverlappedOutputAnno<-cbind(Obj1NDR,rowSums(Obj1Obj2NDRoverlapped))
	Obj2NDRoverlappedOutputAnno<-cbind(Obj2NDR,rowSums(t(Obj1Obj2NDRoverlapped)))
	
	
	
	
	Obj1_Obj2<-Obj1NDRoverlappedOutputAnno[Obj1NDRoverlappedOutputAnno[,15]!=0,]
	
	Obj1_only<-Obj1NDRoverlappedOutputAnno[Obj1NDRoverlappedOutputAnno[,15]==0,]
	Obj2_only<-Obj2NDRoverlappedOutputAnno[Obj2NDRoverlappedOutputAnno[,15]==0,]
	
	
	
	write.table(Obj1_Obj2[,1:3],"Obj1_Obj2.bed",quote=F,sep="\t",row.names=F,col.names=F)
	
	
	write.table(Obj1_only[,1:3],"Obj1_only.bed",quote=F,sep="\t",row.names=F,col.names=F)
	write.table(Obj2_only[,1:3],"Obj2_only.bed",quote=F,sep="\t",row.names=F,col.names=F)
	
	
	
	rangeObj1_Obj2<-BED2RangedData(Obj1_Obj2[,1:3],header=FALSE)
	
	rangeObj1_only<-BED2RangedData(Obj1_only[,1:3],header=FALSE)
	rangeObj2_only<-BED2RangedData(Obj2_only[,1:3],header=FALSE)
	
	
	annObj1_Obj2<-annotatePeakInBatch(rangeObj1_Obj2,AnnotationData=TSS.human.GRCh37, featureType=c("TSS"))
	
	annObj1_only<-annotatePeakInBatch(rangeObj1_only,AnnotationData=TSS.human.GRCh37, featureType=c("TSS"))
	annObj2_only<-annotatePeakInBatch(rangeObj2_only,AnnotationData=TSS.human.GRCh37, featureType=c("TSS"))
	
	
	fullAnnObj1_Obj2<-addGeneIDs(annObj1_Obj2,"org.Hs.eg.db",IDs2Add=c("symbol","genename"),feature_id_type="ensembl_gene_id")
	
	
	fullAnnObj1_only<-addGeneIDs(annObj1_only,"org.Hs.eg.db",IDs2Add=c("symbol","genename"),feature_id_type="ensembl_gene_id")
	fullAnnObj2_only<-addGeneIDs(annObj2_only,"org.Hs.eg.db",IDs2Add=c("symbol","genename"),feature_id_type="ensembl_gene_id")
	
	
	write.table(cbind(Obj1_Obj2[,1:3],annObj1_Obj2$feature),"Obj1_Obj2.annotate.bed",quote=F,sep="\t",row.names=F,col.names=F)
	
	write.table(cbind(Obj1_only[,1:3],annObj1_only$feature),"Obj1_only.annotate.bed",quote=F,sep="\t",row.names=F,col.names=F)
	write.table(cbind(Obj2_only[,1:3],annObj2_only$feature),"Obj2_only.annotate.bed",quote=F,sep="\t",row.names=F,col.names=F)
	
	
	write.table(cbind(Obj1_Obj2[,1:3],fullAnnObj1_Obj2$feature,fullAnnObj1_Obj2$symbol,fullAnnObj1_Obj2$genename),"Obj1_Obj2.annotate.bed",quote=F,sep="\t",row.names=F,col.names=F)
	
	write.table(cbind(Obj1_only[,1:3],fullAnnObj1_only$feature,fullAnnObj1_only$symbol,fullAnnObj1_only$genename),"Obj1_only.annotate.bed",quote=F,sep="\t",row.names=F,col.names=F)
	write.table(cbind(Obj2_only[,1:3],fullAnnObj2_only$feature,fullAnnObj2_only$symbol,fullAnnObj2_only$genename),"Obj2_only.annotate.bed",quote=F,sep="\t",row.names=F,col.names=F)
	
	
	if(!is.null(anno)){
		
		
		Obj1NDRoverlappedOutputAnno<-cbind(Obj1NDRanno,rowSums(Obj1Obj2NDRoverlappedAnno))
		Obj2NDRoverlappedOutputAnno<-cbind(Obj2NDRanno,rowSums(t(Obj1Obj2NDRoverlappedAnno)))
		
#CGI TSS
		
		Obj1_Obj2<-Obj1NDRoverlappedOutputAnno[Obj1NDRoverlappedOutputAnno[,16]!=0,]
		
		Obj1_only<-Obj1NDRoverlappedOutputAnno[Obj1NDRoverlappedOutputAnno[,16]==0,]
		Obj2_only<-Obj2NDRoverlappedOutputAnno[Obj2NDRoverlappedOutputAnno[,16]==0,]
		
		
		write.table(Obj1_Obj2[,1:3],"Obj1_Obj2.CGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
		
		write.table(Obj1_only[,1:3],"Obj1_only.CGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
		write.table(Obj2_only[,1:3],"Obj2_only.CGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
		
		
		
		rangeObj1_Obj2<-BED2RangedData(Obj1_Obj2[,1:3],header=FALSE)
		
		rangeObj1_only<-BED2RangedData(Obj1_only[,1:3],header=FALSE)
		rangeObj2_only<-BED2RangedData(Obj2_only[,1:3],header=FALSE)
		
		
		annObj1_Obj2<-annotatePeakInBatch(rangeObj1_Obj2,AnnotationData=TSS.human.GRCh37, featureType=c("TSS"))
		
		annObj1_only<-annotatePeakInBatch(rangeObj1_only,AnnotationData=TSS.human.GRCh37, featureType=c("TSS"))
		annObj2_only<-annotatePeakInBatch(rangeObj2_only,AnnotationData=TSS.human.GRCh37, featureType=c("TSS"))
		
		
		fullAnnObj1_Obj2<-addGeneIDs(annObj1_Obj2,"org.Hs.eg.db",IDs2Add=c("symbol","genename"),feature_id_type="ensembl_gene_id")
		
		fullAnnObj1_only<-addGeneIDs(annObj1_only,"org.Hs.eg.db",IDs2Add=c("symbol","genename"),feature_id_type="ensembl_gene_id")
		fullAnnObj2_only<-addGeneIDs(annObj2_only,"org.Hs.eg.db",IDs2Add=c("symbol","genename"),feature_id_type="ensembl_gene_id")
		
		
		write.table(cbind(Obj1_Obj2[,1:3],annObj1_Obj2$feature),"Obj1_Obj2.annotate.CGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
		
		write.table(cbind(Obj1_only[,1:3],annObj1_only$feature),"Obj1_only.annotate.CGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
		write.table(cbind(Obj2_only[,1:3],annObj2_only$feature),"Obj2_only.annotate.CGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
		
		
		write.table(cbind(Obj1_Obj2[,1:3],fullAnnObj1_Obj2$feature,fullAnnObj1_Obj2$symbol,fullAnnObj1_Obj2$genename),"Obj1_Obj2.annotate.CGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
		
		write.table(cbind(Obj1_only[,1:3],fullAnnObj1_only$feature,fullAnnObj1_only$symbol,fullAnnObj1_only$genename),"Obj1_only.annotate.CGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
		write.table(cbind(Obj2_only[,1:3],fullAnnObj2_only$feature,fullAnnObj2_only$symbol,fullAnnObj2_only$genename),"Obj2_only.annotate.CGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
		
		if(overlap==0){
			#Non_CGI TSS
			
			Obj1_Obj2<-Obj1NDRoverlappedOutputAnno[Obj1NDRoverlappedOutputAnno[,16]!=0,]
			
			
			Obj1_only<-Obj1NDRoverlappedOutputAnno[Obj1NDRoverlappedOutputAnno[,16]==0,]
			Obj2_only<-Obj2NDRoverlappedOutputAnno[Obj2NDRoverlappedOutputAnno[,16]==0,]
			
			
			write.table(Obj1_Obj2[,1:3],"Obj1_Obj2.nonCGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
			
			write.table(Obj1_only[,1:3],"Obj1_only.nonCGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
			write.table(Obj2_only[,1:3],"Obj2_only.nonCGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
			
			
			rangeObj1_Obj2<-BED2RangedData(Obj1_Obj2[,1:3],header=FALSE)
			
			rangeObj1_only<-BED2RangedData(Obj1_only[,1:3],header=FALSE)
			rangeObj2_only<-BED2RangedData(Obj2_only[,1:3],header=FALSE)
			
			
			annObj1_Obj2<-annotatePeakInBatch(rangeObj1_Obj2,AnnotationData=TSS.human.GRCh37, featureType=c("TSS"))
			
			annObj1_only<-annotatePeakInBatch(rangeObj1_only,AnnotationData=TSS.human.GRCh37, featureType=c("TSS"))
			annObj2_only<-annotatePeakInBatch(rangeObj2_only,AnnotationData=TSS.human.GRCh37, featureType=c("TSS"))
			
			
			fullAnnObj1_Obj2<-addGeneIDs(annObj1_Obj2,"org.Hs.eg.db",IDs2Add=c("symbol","genename"),feature_id_type="ensembl_gene_id")
			
			fullAnnObj1_only<-addGeneIDs(annObj1_only,"org.Hs.eg.db",IDs2Add=c("symbol","genename"),feature_id_type="ensembl_gene_id")
			fullAnnObj2_only<-addGeneIDs(annObj2_only,"org.Hs.eg.db",IDs2Add=c("symbol","genename"),feature_id_type="ensembl_gene_id")
			
			
			write.table(cbind(Obj1_Obj2[,1:3],annObj1_Obj2$feature),"Obj1_Obj2.annotate.nonCGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
			
			write.table(cbind(Obj1_only[,1:3],annObj1_only$feature),"Obj1_only.annotate.nonCGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
			write.table(cbind(Obj2_only[,1:3],annObj2_only$feature),"Obj2_only.annotate.nonCGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
			
			
			write.table(cbind(Obj1_Obj2[,1:3],fullAnnObj1_Obj2$feature,fullAnnObj1_Obj2$symbol,fullAnnObj1_Obj2$genename),"Obj1_Obj2.annotate.nonCGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
			
			write.table(cbind(Obj1_only[,1:3],fullAnnObj1_only$feature,fullAnnObj1_only$symbol,fullAnnObj1_only$genename),"Obj1_only.annotate.nonCGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
			write.table(cbind(Obj2_only[,1:3],fullAnnObj2_only$feature,fullAnnObj2_only$symbol,fullAnnObj2_only$genename),"Obj2_only.annotate.nonCGI.bed",quote=F,sep="\t",row.names=F,col.names=F)
			
		}
		
		
		
	}
}


